############################################
# GHS 
# Date: MAY 19 2014
# This version updates heterogeneous effects analysis in main paper to regression (see changes to table 6 and 7)
# and updates the appendix to include RI approach to hetero effects as a sensitivity analysis to both 
# Price Induced Flattening tests presented in the online appendix
############################################

#############################################
# TABLE OF CONTENTS
#############################################

#   MAIN ANALYSIS
#   Table 2     Access to Existing Communication Channels
#   Table 3     Politically Engaged
#   Figure 1	  Traditional political engagement, wealth and age as a function of measures of access and marginalization	
#   Figure 2    Flattening	
#   Table 4	  	Flattening Participation, test of H1	
#   Figure 3	  Sensitivity of Results 	
#   Figure 4 	  Patterns of Representativeness of Messaging, H1.2		
#   Table 5	  	Price Effects	
#   Table 6   	Price Induced Flattening 1: Test of H 3.1	
#   Table 7	  	Price Induced Flattening 2: Test of H 3.2
#   Table 8	  	Types of messages		
#   Figure 5  	Message Content and Pricing

#######################################################################################################
#                           MAIN ANALYSIS
#######################################################################################################


#######################################################################################################
# Table 2: Access
#######################################################################################################


acc.tbl <- matrix(NA, 16, 6)

vars = cbind(D$q44, D$q47, D$computer, D$travel, D$geo_access, D$ACCESS)
vars2 = cbind(D$POORBI, D$GENDER, D$COG_SMP, D$COE_SMP)


for(z in 1:ncol(vars2)){
  
  for(v in 1:ncol(vars)) {
    tst <- t.test(vars[,v]~vars2[,z], var.equal=FALSE)
    acc.tbl[((z*4)-3):((z*4)-2),v] <- as.matrix(round(tst$estimate, 3))
    acc.tbl[((z*4)-1), v] <- acc.tbl[((z*4)-3), v] - acc.tbl[((z*4)-2), v]
    acc.tbl[(z*4),v] <- round(tst$p.value,3)
  }       
}

colnames(acc.tbl) <- c("Phone Access", "SMS Access", "Computer Access", "Travel Outside Village", "Proximity", "Access Index")
rownames(acc.tbl) <- c("Richer Half", "Poorer Half", "Difference1", "P value1",  "Female", "Male", "Difference2", "P value2",
                       "Noncogender", "Cogender","Difference3", "P value3", "Non-coethnic", "Coethnic", "Difference4", "P value4")
xtable(acc.tbl)


#######################################################################################################
# Table 3: Political Engagement
#######################################################################################################\
eng.tbl <- matrix(NA, 16, 11)

vars = cbind(D$engaged_party, D$engaged_lc1, D$community_meeting, D$raise_issue, D$demonstration,  D$rally,
             D$letter,  D$q56_new,  D$q66bb,  D$q119_SMP,  D$ENGBI)

vars2 = cbind(D$POORBI, D$GENDER, D$COG_SMP, D$COE_SMP)

for(z in 1:ncol(vars2)){
  
  for(v in 1:ncol(vars)) {
    tst <- t.test(vars[,v]~vars2[,z], var.equal=FALSE)
    eng.tbl[((z*4)-3):((z*4)-2),v] <- as.matrix(round(tst$estimate, 3))
    eng.tbl[((z*4)-1), v] <- eng.tbl[((z*4)-3), v] - eng.tbl[((z*4)-2), v]
    eng.tbl[(z*4),v] <- round(tst$p.value,3)
  }       
}

colnames(eng.tbl) <- c("Part Member", "Village Committee", "Community Meetings", "Raise issue", 
                       "Protest", "Election Rally", "Write Letter", "MP Meeting", "Voted MP elections", "Talked to MP", "Engaged Index")
rownames(eng.tbl) <- c( "Richer Half", "Poorer Half", "Difference1", "P value1",  "Female", "Male", "Difference2", "P value2",
                       "Noncogender", "Cogender", "Difference3", "P value3", "Non-coethnic", "Coethnic", "Difference4", "P value4")
xtable(eng.tbl)




#######################################################################################################
# Figure 1: Traditional political engagement wealth and age as a function of measures of access and
# marginalization
#######################################################################################################
pdf(file="../L/figures/measuresNEWBLOCKS.pdf", width=10, height=9) 
	deg=3
	par(mfrow=c(2,2))
	lgraph(D$ACCESS_RANK, 	D$ENGAGED, TRUE, main="Political engagement given access", xlab="Access  (percentile)", ylab="Political Engagement", deg=deg, ylim=c(-.5, .5))
	lgraph(D$MARG_RANK, 	D$ENGAGED, TRUE, main="Political engagement given marginalization", xlab="Marginalization Index (Percentile)", ylab="Political Engagement", deg=deg, ylim=c(-.5, .5))
	lgraph(D$WEALTH_RANK, 	D$ENGAGED, TRUE, main="Political engagement given wealth", xlab="Wealth (percentile)", ylab="Political Engagement", deg=deg, ylim=c(-.5, .5))
	lgraph(D$AGE, D$ENGAGED, D$AGE>=16 & D$AGE<=80, main="Political engagement given  age", xlab="Age", ylab="Political Engagement", deg=deg, ylim=c(-.5, .5))
dev.off()


#######################################################################################################
# Figure 2  Flattening	
#######################################################################################################

pdf(file="../L/figures/marginalizedlfrNEWBLOCKS.pdf", width=6, height=5) 
	deg=3
	lgraph(D$MARG_RANK,	D$SMS,	D$V6_SMS==1, main="ICT engagement given marginalization", xlab="Marginalization Index (Percentile)", ylab="SMS takeup", deg=deg, ylim=c(0, .1))
dev.off()



#######################################################################################################
# Table 4	  Flattening Participation, test of H1	
#######################################################################################################

# Main Table: H1:
xtable(h1tabf(D$V6_SMS==1, Y=D$MARGINAL, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share Marginalized"))


	# Manual check
	.pctile =  1- mean(DS$SMS)				# Cutoff
	.ENGAGED = (rank(DS$ENGAGED)> max(rank(DS$ENGAGED))*.pctile) 
	mean((DS$MARGINAL)[.ENGAGED])
	mean((DS$MARGINAL)[DS$SMS==1])
	mean((DS$MARGINAL)[DS$SMS==1 & DS$PRICE==0])
	mean((DS$MARGINAL)[DS$SMS==1 & DS$PRICE==1])

	.pctile2 =  1- mean(DS$SMS[DS$PRICE==2])				# Cutoff
	.ENGAGED2 = (rank(DS$ENGAGED)> max(rank(DS$ENGAGED))*.pctile2) 
	mean((DS$MARGINAL)[.ENGAGED2])
	mean((DS$MARGINAL)[.ENGAGED2 & DS$PRICE==2])  # Note that among those in the high price group there are more marginalized subjects; by chance
	mean((DS$MARGINAL)[DS$SMS==1 & DS$PRICE==2])



#######################################################################################################
# Figure 3   Sensitivity of H1 Results
#######################################################################################################

  h1tabf2   = function(C, Y = D$MARGINAL, YLAB = "Share Marginalized", dig=3, PCT){
    ENGB 	= (rank(D$ENGAGED[C])>max(rank(D$ENGAGED[C]))*PCT)   # DEFINE ENG VARIABLE FOR THE SET IN QUESTION
    SMSC  = D$SMS[C]
    msur 	<- systemfit(list(SMSreg = Y[C]~SMSC, ENGreg = Y[C]~ENGB), data=D)
    CM	=  coef(summary(msur))
    LH0   = linearHypothesis(msur, "SMSreg_SMSC + SMSreg_(Intercept) - ENGreg_ENGBTRUE - ENGreg_(Intercept)", test = "Chisq", na.rm=TRUE)
    X= cbind(round(CM[1,1]+CM[2,1]- CM[3,1]- CM[4,1], dig), round(LH0[[4]][[2]], dig))
    return(X)
  }

# Uptake rate by definition of "engaged"


  h1tabf3   = function(PCT){
    ENGB 	= (rank(DS$ENGAGED)>max(rank(DS$ENGAGED))*PCT)   # DEFINE ENG VARIABLE FOR THE SET IN QUESTION
    log(mean(ENGB[DS$MARGINAL==0]) / mean(ENGB[DS$MARGINAL==1])) 
	}

  h1tabf3(.5)
  h1tabf3(.8)
  
  
grain = 200
popsplits <- seq(.02, .98, length=grain)
results <- matrix(NA, grain, 2)
for(i in 1:grain){results[i,] <- h1tabf2(D$V6_SMS==1, PCT=popsplits[i])}

pdf(file="../L/figures/FlatteningSensitivityNEWBLOCKS.pdf", height=6.5, width=7.5)
	
plot(0,0, xlim= c(0,1), ylim=c(-.1,.4), xlab="Threshold for classification as `engaged'", ylab="Estimated Flattening", type="n")
	points(popsplits, results[,2], type="l", col="grey")
		text(.08, 0.24, ("p-value"), col="grey") 
	abline(0,0)
	abline(h=.05, col= "red")
	abline(h=.1, col= "red")
	abline(v=.955, col= "black", lty=2)
	points(popsplits, results[,1], type="l", col="black", lwd=1.5)
	axis(4, at=c(.05, .1),
     col.axis="red", las=2, cex.axis=0.7, tck=-.01)


dev.off()


#######################################################################################################
# Figure 4   Patterns of Representativeness of Messaging, H1.2		
#######################################################################################################

A = table(DS$priority)/sum(table(DS$priority))
E = table(DS$priority[DS$ENG==1])/sum(table(DS$priority[DS$ENG==1]))
S = table(DS$priority[DS$SMS==1])/sum(table(DS$priority[DS$SMS==1]))
N = as.table(cbind(length(DS$priority),length(DS$priority[DS$ENG==1]),length(DS$priority[DS$SMS==1]) ))
Z=rbind(A, E)
T=rbind(A, S)

# TABLE OF MESSAGES  

ISSUET = round(cbind(A, E, S),3)
ISSUET = rbind(ISSUET, N)
rownames(ISSUET) <- c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other", "Observations")
colnames(ISSUET) <- c("All subjects", "Engaged subjects", "SMS senders")
xtable(ISSUET)

(A-E)%*%(A-E)/2
(A-S)%*%(A-S)/2

def.par <- par(no.readonly = TRUE)
par(def.par)

pdf(file="../L/figures/PriorityDistNEWBLOCKS.pdf", width=13, height=8)
	par(mfrow=c(1,2))
	par(xaxt="n")
	colors=c("grey85","grey20")
	enghist <- barplot(Z[1:2,], beside = TRUE,
					   ylim = c(0, max(T)), yla="Group Share", border="blue")
	legend("topright", c("Engaged", "Population"), fill=colors, cex=.95)
	lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
	axis(1, at=seq(1, 39, by=3), labels = FALSE)
	text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)

	colors=c("grey55","grey20")
	smshist <- barplot(T[1:2,], beside = TRUE,
					   ylim = c(0, max(T)), yla="Group Share", border="blue")
	legend("topright", c("SMS", "Population"), fill=colors, cex=.95)
	lablist=as.vector(c("Welfare", "Economy", "Taxes", "Credit", "Infrastructure", "Agriculture", "Education", "Health", "Water", "Electricity", "Crime", "Governance", "Other"))
	axis(1, at=seq(1, 39, by=3), labels = FALSE)
	text(seq(1, 39, by=3), par("usr")[3]-0.015, labels = lablist, srt = 90, pos = 1, xpd = TRUE, cex=.95, font=2)
	dev.off()
par(def.par)


#######################################################################################################
# Table 5	  Price Effects	
#######################################################################################################

# Basic price effect; lm version (reports 2 sided)
# Note higher "price" means more subsidization
# Note that I changed the code a little here since Y_REDU was not defiend
Y_REDU=cbind(DS$SMS, DS$SMS_B, DS$SMS_V)

for(y in 1:3){
  print(coef(summary(lm(Y_REDU[,y]~DS$PRICE+as.factor(DS$BLOCK))))[2,], digits=3)
}

#top part of the table
MainT=cbind(
  sapply(0:2, function(i) mean(DS$SMS[DS$PRICE==i])),
  sapply(0:2, function(i) mean(DS$SMS_B[DS$PRICE==i])),
  sapply(0:2, function(i) mean(DS$SMS_V[DS$PRICE==i])),
  sapply(0:2, function(i) sum(DS$PRICE==i, na.rm=T)))

rownames(MainT) <- c("Full", "Subsidy", "Free")
colnames(MainT) <- c("Any SMS", "Public SMS", "Private SMS", "Obs")
xtable(MainT, digits=3)

# BOTTOM PART OF THE TABLE
# DEFINE REDUCED Y
Y_REDU=cbind(DS$SMS, DS$SMS_B, DS$SMS_V)
OUT = matrix(NA, 12,3)
for(y in 1:3){
  OUT[1:3,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==1), DS$BLOCK, PM, C = (DS$PRICE!=2), addn=TRUE, round=3, onesided=3)  # T_FS  
  OUT[4:6,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=0), addn=TRUE, round=3, onesided=3)  # T_SO
  OUT[7:9,y]     <-ri_lm(Y_REDU[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=1), addn=TRUE, round=3, onesided=3)  # T_FO
}
for(y in 1:3){
  OUT[10:12,y]   <-ri_lm(Y_REDU[,y],  DS$PRICE,  DS$BLOCK, PM, addn=TRUE, round=3, onesided=3)   # Linear   
}

# H4 test using SUR
msurH4 	<- systemfit(list(MV = SMS_V~PRICE + as.factor(BLOCK), MB = SMS_B~PRICE + as.factor(BLOCK)), data=DS)
CM_H4	=  coef(summary(msurH4))

###coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"]

LH_H4 	= linearHypothesis(msurH4, "MB_PRICE  - MV_PRICE", test = "Chisq", na.rm=TRUE)
pH4_plus = round(LH_H4[[4]][2]/2,3)
pH4_plus

OUTp= rbind(cbind(OUT, c(rep("",11),pH4_plus)), c("H2 test","","",""))
C1 = c("Subsidy vs. Full Price", "", "","Free vs. Subsidy","", "","Free vs. Full Price","", "", "Linear Trend", "","", "")
C2 = c(rep(c("ATE", paste("(p)"), "(N)"),3), c("Trend", "(p)", "(N)"), "")
OUT2 = cbind(C1,C2, OUTp)
OUT2[11,6]<-round(coef(msurH4$eq[[1]])["PRICE"]- coef(msurH4$eq[[2]])["PRICE"],3)

colnames(OUT2)<-c("Treatment", "Effect", "Any", "Public","Private", "H4 test")
print(xtable(OUT2, digits=3), include.rownames=FALSE)


#######################################################################################################
# Table 6 	Price Induced Flattening 1: Test of H 3.1	
#######################################################################################################
OUT6 <- matrix(0, 6, 4)
DS$RICHBI = ifelse(DS$POORBI==0, 1,0)
DS$NACC = ifelse(DS$ACC==0, 1,0)

#diff in marginal effects
poverty.diff  <- lm(SMS ~ (POORBI*PRICE) + PRICE + POORBI + ACC + as.factor(BLOCK), data=DS)
poverty.diff2  <- lm(SMS ~ (RICHBI*PRICE) + PRICE + RICHBI + ACC + as.factor(BLOCK), data=DS)

s=summary(poverty.diff) 
OUT6[2,1]     <- coef(poverty.diff)["PRICE"]
OUT6[2,2]     <- coef(poverty.diff2)["PRICE"]
OUT6[2,4]     <- coef(poverty.diff)["POORBI:PRICE"] 
OUT6[3,1]     <- summary(poverty.diff)$coefficients["PRICE",2]
OUT6[3,2]     <- summary(poverty.diff2)$coefficients["PRICE",2]
OUT6[3,4]     <- pt(coef(s)[, 3], s$df[2], lower=FALSE)["POORBI:PRICE"] 

#diff in marginal effects
acc.diff  <- lm(SMS ~ (ACC*PRICE) + PRICE + ACC + POORBI + as.factor(BLOCK), data=DS)
acc.diff2  <- lm(SMS ~ (NACC*PRICE) + PRICE + NACC + POORBI + as.factor(BLOCK), data=DS)

m=summary(acc.diff)
OUT6[5,1] <- coef(acc.diff)["PRICE"]
OUT6[5,2] <- coef(acc.diff2)["PRICE"]
OUT6[5,4] <- coef(acc.diff)["ACC:PRICE"]
OUT6[6,1]     <- summary(acc.diff)$coefficients["PRICE",2]
OUT6[6,2]     <- summary(acc.diff2)$coefficients["PRICE",2]
OUT6[6,4] <- pt(coef(m)[, 3], m$df[2], lower=FALSE)["ACC:PRICE"]

OUT6 <- round(OUT6, 3)

rows = c("", "Marginal Effect of Subsidy by Poverty","","", "Marginal Effect of Subsidy by Political Access", "")
OUT6 <- cbind(rows, OUT6)
OUT6[1, 1:5] <- c("", "Rich", "Poor", "Hypothesized Difference", "Above-Below")
OUT6[4, 1:5] <- c("", "Low Access", "High Access", "Hypothesized Difference", "Above-Below")

xtable(OUT6)

#######################################################################################################
# Table 7	  Price Induced Flattening 2: Test of H 3.2	
#######################################################################################################

# Difference Marginalized v Non-Marginalized, no controls
DS$NONMARGINAL = ifelse(DS$MARGINAL==0, 1,0)
marg.diff  <- lm(SMS ~ (PRICE*MARGINAL) 	+ PRICE + MARGINAL + as.factor(BLOCK), data=DS)
marg.diff2 <- lm(SMS ~ (PRICE*NONMARGINAL) 	+ PRICE + NONMARGINAL + as.factor(BLOCK), data=DS)

v=summary(marg.diff)
q=summary(marg.diff2)

OUT7 <- matrix(NA, 3 ,2)
	OUT7[1,1] <- coef(marg.diff2)["PRICE"]
	OUT7[1,2] <- pt(coef(q)[, 3], q$df[2], lower=FALSE)["PRICE"]
	OUT7[2,1] <- coef(marg.diff)["PRICE"]
	OUT7[2,2] <- pt(coef(v)[, 3], v$df[2], lower=FALSE)["PRICE"]

	OUT7[3,1] <- coef(marg.diff)["PRICE:MARGINAL"]
	OUT7[3,2] <- pt(coef(v)[, 3], v$df[2], lower=FALSE)["PRICE:MARGINAL"]
	OUT7 <- round(OUT7, 3)

	OUT7 <- cbind(c("Marginal Effect of Price on Marginalized", "Marginal Effect of Price on Non-marginalized", "Difference"), OUT7)
xtable(OUT7)

SEND = DS$SMS==1
h(DS$MARGINAL[SEND], DS$PRICE[SEND])
bss=bsh(DS$MARGINAL[SEND],PM[SEND,])
hist(bss)

# changed to reflect treatment as subsidy (0 is full price, 1 is subsidy, 2 is complete subsidy)
h5.3tab = function(Y, X, name, PERM=PM[SEND,]) (c(
  paste("Share of ",  name, " among high price senders=", round(mean(Y[X==0]),3), sep=""), 
  paste("n_0  =", round(length(Y[X==0]),3)),
  paste("Share of ",  name, " among mid price senders= ",round(mean(Y[X==1]),3), sep=""),
  paste("n_S  =", round(length(Y[X==1]),3)),
  paste("Share of ",  name, " among low price senders= ",round(mean(Y[X==2]),3), sep=""),
  paste("n_F  =", round(length(Y[X==2]),3)),
  paste("Trend from high to low=",  round(h(Y, X),3)), 
  paste("(p (neg))=",mean(h(Y, X)>=bsh(Y, PERM))),
  paste("(p (two sided))=",mean(abs(h(Y, X))<=abs(bsh(Y, PERM)))),
  paste("(p (pos))=",mean(h(Y, X)<=bsh(Y, PERM)))
))
h5.3tab(DS$MARGINAL[SEND],DS$PRICE[SEND], "marginal respondents")
#h5.3tab(DS$POORBI[SEND],DS$PRICE[SEND], "poor respondents")
#h5.3tab(DS$WOMEN[SEND],DS$PRICE[SEND], "women")

#######################################################################################################
# Table 8	  Types of messages		
#######################################################################################################
table(DS$PubPriv, DS$SMS)  

#######################################################################################################
# Figure 5	Message Content and Pricing		
#######################################################################################################
#treatment as subsidy (0 is full price, 1 is subsidy, 2 is complete subsidy)

counts <- as.matrix(table(DS$PRICE, DS$PubPriv))
counts <-counts[,-1] #remove uncoded messages

fullprice <- sum(DS$SMS[DS$PRICE==0 & DS$PubPriv!=0], na.rm=TRUE)
counts1 <- counts[1,]/fullprice
halfprice <- sum(DS$SMS[DS$PRICE==1 & DS$PubPriv!=0], na.rm=TRUE)
counts2 <- counts[2,]/halfprice
freeprice <- sum(DS$SMS[DS$PRICE==2 & DS$PubPriv!=0], na.rm=TRUE)
counts3 <- counts[3,]/freeprice

counts <- rbind(counts1,counts2,counts3)

pdf(file="../L/figures/PriceEffects2NEWBLOCKS.pdf", width=7, height=5.5) 
label <- c("Full Price", "Partially Subsidized", "Free")
plot <- barplot(counts, main="Messaging by Type and Price",
        ylim=c(0, .7),
        xlab="Type of Message", ylab="Proportion of messages within each price group",
        col=c("lightgrey", "darkgrey", "black"),
        legend = label,
        names.arg = c("Personal", "Local", "Constituency", "Group", "National"),
        cex.names=1,
        beside=TRUE)
dev.off()


